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1.  Introduction 


A  physical  restriction  in  most  conventional  continuum  models  is  their  inability  of  taking  into 
account  sufficient  microscopic  degrees  of  freedom  in  the  constitutive  framework.  For  crystalline 
metallic  materials,  the  dislocation-related  microstructural  degrees  of  freedom  responsible  for 
dislocation  substructure  evolution  are  critically  important  to  model  and  simulate  multiscale 
aspects  of  plastically  deforming  metals.  A  solution  to  this  problem  has  been  recently  provided  by 
the  Field  Theory  of  Multiscale  Plasticity  (FTMP)  advocated  by  Hasebe  {1-4).  In  this  relatively 
new  theoretical  framework,  additional  degrees  of  freedom  are  implemented  via  FTMP-based 
tensor  of  incompatibility  derived  from  the  differential  geometric  curvature  tensor  constructed 
from  gradients  of  the  metric  tensor  of  the  plastically  deformed  (or  elastically  unloaded) 
incompatible/anholonomic  crystalline  space  (4-6),  to  be  ultimately  incorporated  in  a  crystalline 
plasticity-based  constitutive  framework. 

In  this  work,  we  propose  a  new  model  for  deformation  twinning  based  on  FTMP.  Twinning  is  an 
alternative  mode  of  plastic  deformation  to  slip  via  dislocation  glide,  the  former  involving 
collective  motion  of  partial  dislocations.  Deformation  twinning  is  a  major  deformation  mode  for 
hexagonal  close-packed  (HCP)  metals  (here  HCP  refers  to  any  hexagonal  metal,  not  necessarily 
one  with  the  ideal  c/a  ratio  for  closest  packing,  where  c  and  a  are  lattice  parameters  of  the 
hexagonal  unit  cell)  and  also  for  some  face-centered  cubic  (FCC)  metals  with  low  stacking  fault 
energy  (SFE)  (e.g.,  less  than  25  mJ/m2)  and  in  body-centered  cubic  (BCC)  metals  under  high 
strain  rates  and/or  low  temperatures.  Even  in  FCC  metals  with  intermediate  SFE  around 
50-70  mJ/m2,  such  as  copper  (Cu)  and  nickel  (Ni),  where  most  plastic  deformation  can  be  carried 
by  glide  of  dislocations,  mechanical  twinning  can  still  take  place,  e.g.,  under  hypervelocity  impact 
loading.  Since  both  twinning  and  slip  resistance  can  be  affected  by  densities  of  continuously 
distributed  dislocations  introduced  via  accommodation  of  excessive  local  deformation,  extended 
use  of  the  incompatibility  tensor-based  model  is  expected  be  effective  for  their  descriptions. 

In  this  report,  a  brief  theoretical  overview  and  a  “flow-evolutionary  working  hypothesis”  are 
given  in  section  2.  The  constitutive  framework  is  outlined  in  sections  3  (slip)  and  4  (twinning). 
Representative  results  from  deformation  simulations  of  cubic  and  hexagonal  metallic  crystals 
(Cu  and  magnesium  [Mg],  respectively)  are  presented  in  sections  5  and  6.  Bold  font  is  used  for 
vectors  and  tensors  and  italic  for  scalars  with  summation  over  repeated  indices. 
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2.  Theoretical  Background 


2.1  Differential  Geometric  Description  of  Defect  Fields 

Let  V  denote  the  covariant  derivative  with  respect  to  arbitrary  connection  r .  In  spatial 
coordinates  x'  (i  =  1,2,3),  the  covariant  derivative  of  a  smooth  vector  field  F(x)  is 

vyj  =dlvJ+r/kvk.  (i) 


The  covariant  derivative  of  the  space  itself  can  be  defined  by 

(Vx)'  =  (dkx‘  +  n.x1  )dxl  =  dx1  +  n.x1  dxk . 


Let  indices  in  brackets  be  skew,  indices  in  parentheses  be  symmetric,  and  indices  in  vertical  bars 
be  excluded  from  such  operations,  e.g., 

^  Ay]  =  A/  -  A'i  ’  U)  =  Ay  +  A i  ’  2BW\k]  =  Byk  ~  & kji  •  (^) 

Torsion  and  curvature  tensors  are  defined,  respectively,  for  a  generic  connection  as 

S’kl=r’uy  x„„"=2(a, (4) 


[k-1  l]m  1  [Jfc|p|J  7], 


When  torsion  S  is  nonzero,  coefficients  of  /^describe  a  non-Riemannian  space.  When  curvature 
R  vanishes,  the  space  is  said  to  be  flat. 

In  finite  deformation  plasticity,  a  multiplicative  decomposition  of  the  deformation  gradient  F 
into  elastic  and  plastic  terms  is  standard: 

F  =dx/dX=Fe -Fp  ^l  +  pe  +pp ,  (5) 

where  in  the  final  approximation,  small  deformations  are  assumed  such  that,  for  the  displacement 
field  u,  in  Cartesian  coordinates,  elastic  and  plastic  distortions  are  additive: 

duj  /  dXi  =4  duj  /dxi  =  dgUj  =  (3en  +  •  (6) 


Two  connections  are  introduced: 


The  first  has  vanishing  curvature  but  nonvanishing  torsion  S  ,  and  the  second  has  vanishing 
torsion  but  nonvanishing  curvature  R  .  Contractions  of  these  higher-order  tensors  considering 
symmetry  results  in  well-known  second-order  tensors.  They  are  respectively  called  “dislocation 
density  tensor”  and  “incompatibility  tensor,”  given  in  the  geometrically  linear  regime  by  the  curl 
of  the  elastic  or  plastic  distortion  tensor  and  double  curl  of  the  elastic  or  plastic  strain  tensor: 


2 


(8) 


a,]=-smSkj  ^s"ddkP 


ej 

■l 


-sMdkP 


Pj 

•/’ 


^■1  —  1  c?ikl  (?  j,nn  J?  _  pjmnps  ps  „e  _  „ikl  „  jmn ps  ps  „p 

rj  -  4 1,  t,  i\ H  —  t,  t,  ukumt,ln  -t.  t,  ukumt,ln , 


(9) 


where  elastic  and  plastic  strains  are  £eln  =  and  £j’n  =  /3[ln) ,  and  where  £,u  =  g~V2eM  with 
g  =  det[g(/]  =  det[3(x  -o^x]  and  eM  is  the  dimensionless  alternator  symbol.  In  the  right  two 
equalities  of  equation  9,  terms  on  the  order  of  the  square  of  the  connection  coefficients  are 
assumed  small  relative  to  strain  gradients.  Complete  derivations  can  be  found  elsewhere  ( 4 ,  6,  7). 
Henceforth,  Cartesian  coordinates  are  used  ( glj  =  Stj),  and  linearized  tensors  in  equations  8  and  9 
are 


=  ZmdiPkj .  %  =  emejnmdkdm£[n .  (10) 

The  important  point  to  note  here  is  that  these  quantities  are  expressed  as  gradients  of  plastic 
distortion  or  plastic  strain,  meaning  the  theory  intrinsically  requires  strain  gradients  at  least  up  to 
the  second  order. 

2.2  Flow-Evolutionary  Hypothesis  in  FTMP 

Hasebe’s  “flow-evolutionary  hypothesis”  asserts  laws  for  evolution  of  inhomogeneous  fields  and 
attendant  local  plastic  flow  accompanied  by  energy  dissipation.  The  notion  of  “duality”  between 
fluctuating  hydrostatic  stress  and  deviatoric  strain  fields  introduced  in  FTMP  {2-4)  can  be 
embodied  by  this  law,  although  it  still  is  a  “working  hypothesis”  deserving  further  verification 
and  validation.  Specifically,  it  represents  an  interrelationship  between  locally  stored  strain 
energy  (fluctuation  part)  and  local  plastic  flow  (in  the  form  of  incompatibility  tensor  field  77 )  as 

has  been  discussed  in  the  context  of  polycrystalline  plasticity  (8).  A  brief  motivation  for  this 
hypothesis  follows. 

Consider  an  elastic  medium  with  internal  energy  density  W  and  kinetic  energy  density  JC : 

W=W(€,x),  K  =  \pv-v,  C  =  JC-W,  n  =  K+W.  (11) 

Here,  C  and  H  are  the  Lagrangian  and  Hamiltonian  energy  densities,  respectively.  Let  the 
Cauchy  stress  be  a  =  SW/ds  and  the  energy  momentum  tensor  be  (9) 

Tij  =  +  °ikdjuk  X  djTy  +  pd{VjdjU t  )/dt  =  -ft  =  (< 3C/dxt , )  lexp .  (12) 

Here, /is  the  local  force  per  unit  volume  acting  on  an  inhomogeneity  or  defect. 

In  the  static  case,  and  when  W=  W(e) ,  i.e.,  the  body  contains  no  heterogeneities,  then  T  is 
divergence-free.  The  incompatibility  tensor  is  symmetric  and  divergence-free: 

dflij  =  dflji  =  ejkieimP  =  e[jk]ieimn^(Pk)^m£ln  =  (13) 
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Under  these  assumptions,  since  both  energy-momentum  tensor  and  incompatibility  tensor  are 
divergence-free,  the  following  linear  relationship  can  be  postulated: 

Vij  =  Kij kl^kl  ’  Kijkl  =  2  K^ik^ jl  +  &i l^jk )  ’  (^4) 

where  isotropy  is  assumed  in  the  definition  of  the  constants  Kijkl ,  so  that  the  first  of  equation  14 
is 


%  =  Ksrm  =>  a  flu  =  K0,srm  =  o ,  (is> 

consistent  with  the  vanishing  divergence  requirement  if  the  skew  part  of  T  is  divergence-free 
and/or  vanishing.  Equation  15  can  be  inverted  to  give 

srm=fd%,  (16) 

where  //is  the  shear  modulus  and  /  a  length  constant.  Physically,  equation  14  appears  sensible 
since  the  energy-momentum  tensor  is  a  driving  force  for  the  evolution  or  motion  of  defects, 
which,  in  turn,  are  represented  by  incompatibility.  A  complementary  hypothesis  is  (4) 

tri/  cc  SH  =S(1C+ W) ,  (17) 

which  can  be  formulated  using  the  4-D  permutation  tensor  entering  the  generalized  cross  product 
(10).  The  proposed  relationship  between  incompatibility  and  elastic  energy  will  be  further 
explored  in  the  context  of  numerical  results  in  section  5.2. 


3.  Constitutive  Framework 


3.1  Shear  Stress-Shear  Strain  Relationship 

As  a  tentative  vehicle  for  accommodating  the  field  theoretical  concepts  overviewed  in  the 
previous  section  to  practical  applications,  the  kinematical  foundations  of  crystal  plasticity  are 
utilized.  The  constitutive  equation  to  be  used,  on  the  other  hand,  puts  its  basis  on  the  statistical 
mechanics-based  dislocation  dynamics  with  which  wide  ranges  of  strain  rate  and  temperature  can 
be  taken  into  account  in  a  unified  manner. 

The  current  field  theoretical  notions  and  quantities,  including  the  interaction  field  framework, 
can  be  easily  implemented  into  the  crystalline  plasticity  model  through  strain  gradient  terms  ( 8 , 
11, 12).  Hasebe  proposed  a  constitutive  model  applicable  to  FCC,  BCC,  and  HCP  metals  based 
on  statistical  mechanics  and  dislocation  dynamics.  The  rate  of  plastic  distortion  is  the  usual  form 
from  crystal  plasticity  theory: 

Fp  ^jlp  =YJNa_jasa®ma ,  (18) 
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with  sa  and  m"  the  slip  direction  and  plane  normal  for  system  a.  The  explicit  form  for  slip  rates 
due  to  resolved  shear  stress  ra  =  a :  sa®ma  is  ( 8 , 11) 


ya=Asz'a  <1  r'a\  Bs  exp 


-i 


l-(r,a/Ka)P  +C5  ,  r'a={ra -TaP)-na . 


(19) 


Here,  As=  pmbLu,  Bs= AG0  /kBT,  Cs=  BLo/b ,  where  K  and  Q  are  drag  and  back  stress, 
respectively,  responsible  for  isotropic  and  kinematic  hardening;  pm,  L,  v,  b,  and  B  are  mobile 


dislocation  density,  mean  flying  distance  of  dislocations,  modified  Debye  frequency,  the 
magnitude  of  Burgers  vector  and  the  damping  coefficient  due  to,  e.g.,  phonon  drag,  respectively; 
and  A G0  is  activation  energy  for  dislocation  processes  at  T  =  0  K.  Exponents  p  and  q  in  equation 
19  are  parameters  specifying  thermal  obstacles  of  interest,  provided  0  <  p  <  1  and  1  <  q  <  2.  In 
the  standard  case,  p  =  q  =  §  is  used  for  representing  dislocation  processes.  Furthermore, 

(•)  =  i(|.|+.)  is  the  Macaulay  parentheses,  with  expressing  the  effective  stress  for  overcoming 


the  Peierls  process  given  by 


r„  = 


1- 


v^o>3  rj 


“I  VPp 


(20) 


where  gP,y0,  pP, and  qp  are  parameters  for  thermal  activation  via  the  Peierls  overcoming  event. 

This  constitutive  model  can  express  stress-strain  responses  for  FCC,  HCP,  and  BCC  metals  over 
a  wide  range  of  strain  rates  and  temperatures,  including  impact  loading  conditions.  For  BCC, 
typically  Cs-0 ;  for  FCC  and  HCP,  typically  r“  /  Ka—  0 . 


3.2  Drag  Stress  Model 


The  hardening  evolution  models  are  affected  through  drag  stress  K° .  The  interactions  of 
dislocation  against  forests  (or  equivalently,  interactions  among  dislocations  mostly  belong  to 
different  slip  systems)  are  regarded  as  being  primarily  responsible  for  the  additional  hardening. 
The  quality  of  the  additional  hardening,  i.e.,  “history,”  is  characterized  by  the  effective  cell  size 
to  be  defined  as  follows  via  hardening  ratio. 


The  time  evolution  of  the  drag  stress  defines  the  instantaneous  hardening  modulus  and  is 
assumed  to  be  expressed  as  the  following  form: 


^ drref  df 


rJ_\z,P 


I, 


(21) 


where  H(/)  =  drref  /dy  is  the  referential  hardening  modulus  with  no  interaction  among  slip 
systems,  while  Qa/j  denotes  the  hardening  ratio  that  evolves  with  histories  and  interactions.  In 

order  to  express  the  nonlocal  actions  associated  with  the  evolution  of  inhomogeneities  based  on 
dislocation  density  and  incompatibility  fields,  strain  gradient  terms  are  additively  introduced  into 
Qap  as  discussed  in  section  3.3. 
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Dimensionless  hardening  ratio  Qap  expresses  the  ratio  of  the  flow  stress  increase  to  the 

referential  level,  measuring  the  additional  hardening,  and  is  further  used  to  evaluate  the  effective 
cell  size,  because  it,  in  principle,  contains  all  the  information  about  the  dislocation  cells  evolved 
(responsible  for  the  additional  hardening).  Therefore,  we  can  assume  that  the  quantity 
characterizes  an  effective  size  of  dislocation  cells,  defined  here  as  (72) 


d=  c/A 


-1/2 


(22) 


where  d0  represents  the  initial  cell  size,  normally  a  fraction  of  the  initial  grain  size. 


3.3  Hardening  Law  and  Field  Theoretical  Strain  Gradient  Terms 

The  model  incorporates  contributions  of  a  and  r}  through  hardening  ratio  Qufj  as  (77) 

+  F‘  (a, ,)] ,  (23) 

where  faK  represents  the  dislocation  interaction  matrix  and  SK/i  expresses  a  history  matrix 

further  given  as  an  increasing  function  of  plastic  work  done  by  the  effective  stress  that  is 
responsible  for  dislocation  processes: 

S^tanhOV'/W^,),  (nosum).  (24) 

Here,  Wj!  defines  the  plastic  work  done  by  the  effective  stress  for  only  the  dislocation  processes: 

W/  =  J Wfdt  =  |  (tp  -  rpP  )  fdt .  (25) 

Function  Fp  {a13 ',r/p)  =  Ff  (aft)  +  F2  (rf)  expresses  the  field  theoretical  “strain  gradient  term” 
contributions  ( 8 ): 

7 ?fi(afi)  =  (k/pfi)(\afi  I  /b)1'2,  F2^rJ/3)  =  (k/qp)(ld\^  \/bf 2  sgitf),  (26) 

where  pp,  qfi  are  coefficients  related  to  contributions  of  dislocations  and  incompatibility  to  the 
change  in  the  effective  cell  size  d  ,  while  ld  represents  the  characteristic  length  of  the  defect 
field  considered,  e.g.,  ld  =  b  for  dislocation  dipoles  and  ld»  10”6m  for  substructures  like  cells. 
Here,  a,s  and  r/’  are  resolved  components  of  ajJ  and  rj^ ,  respectively: 

if  =[(sfi  xmfi)®sP  +sp  <8>sfi  +sp  (27) 

The  first  and  second  terms  in  brackets  respectively  correspond  to  edge  and  screw  components  of 
dislocation  density.  A  similar  assertion  can  be  made  for  incompatibility,  where  the  third  will  be 
called  disclination  type,  though  not  corresponding  to  disclination  theory  in  (6). 
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4.  Modeling  Deformation  Twinning 


The  elastic -plastic  decomposition  (equation  5)  is  extended  to  allow  for  twinning  deformation  (13): 

F=FFFp^  l  +  fie+fi'+fip.  (28) 

Deformation  for  twinning  is  represented  by  F‘  with  approximate  time  rate 

F1  -  P’  =  Y,Mx_jxsx®mx ,  (29) 

with  twinning  direction  sz ,  habit  plane  normal  mx ,  and  shearing  rate  yx  ,  which  can  be 
expressed  in  terms  of  the  change  in  mass  fraction  of  the  twinned  volume  and  the  quantized 
maximum  twinning  shear  (6, 13).  The  geometrically  nonlinear  version  of  the  model  implements 
the  Jaumann  rate  of  the  Kirchhoff  stress  tensor  r  =  (det  F)a .  In  the  geometrically  linear  regime, 

Cauchy  stress  is  given  by 

^=q*/£=q»(5, cm=e2w/ds-ds‘M.  (30) 

The  strain  rate  due  to  deformation  twinning  is  assumed  to  be  driven  by  the  evolution  of  the 
incompatibility  tensor  field  for  the  twinning  mode,  r/z  =  [(sz  xmz) 0sz] :  jj  ,  i.e., 

r=Z AfVU-  e><?rfF3(i-F3/F3“)-  od 

Here,  F3 {if )  denotes  the  FTMP-based  incompatibility  term  defined  as 

F3  =(V qx){riXl<, /b)m  sgn(^) •  (32) 


5.  Numerical  Simulations:  FCC  Copper  (Cu) 


5.1  Model,  Results,  and  Discussion 

The  FTMP-based  incompatibility  model  corresponding  to  slip  and  twinning  degrees  of  freedom 
has  been  introduced  into  the  hardening  law  in  a  crystalline  plasticity-based  constitutive  equation 
applied  to  FCC  Cu.  The  model  has  been  implemented  by  Hasebe  in  an  original  finite-element 
method  (FEM)  code  for  static  implicit  analysis.  Tension/compression  analyses  are  conducted  on 
a  single  slip-oriented  single  crystal  to  check  basic  capabilities  of  the  model  for  describing 
several  typical  features  of  the  deformation  twinning:  nucleation,  growth  (into,  e.g.,  lenticular 
shapes),  lattice  rotation  (satisfying  the  mirror  symmetry),  the  attendant  stress  responses  to  be 
accompanied  by  softening  with  serrations,  and  energy  redistribution  into  the  twin  degrees  of 
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freedom.  Similar  analyses  are  carried  out  on  multigrain  samples  (composed  of  seven  grains)  and 
larger  polycrystalline  samples  to  systematically  investigate  the  effects  of  geometrically  induced 
inhomogeneity  on  nucleation  and  growth  of  multiple  twins  and  the  attendant  stress-strain 
responses. 

Figure  1  shows  the  domain  analyzed  in  the  present  series  of  finite-element  (FE)  simulations,  with 
24  x  24  cross  triangle  elements  for  plane  stress  conditions.  The  sample  is  pulled  in  a  single  slip 
orientation  [125]  to  up  to  e=  0.20  (i.e.,  20%  tensile  elongation). 


Figure  2  displays  representative  results  exhibiting  nucleation  and  the  following  growth  of  twins 
under  tension,  together  with  the  stress  response  accompanied  by  serration.  The  twinned  regions 
yield  lattice  rotations  that  tend  to  saturate  at  around  20°.  Contours  denote  twinned  regions  by 
frequency  of  meeting  a  critical  resolved  shear  stress  condition. 

The  twinning  model  becomes  operative  when  the  critical  condition  for  the  resolved  shear  stress 
is  met,  i.e.,  tz  >  r'cr .  Since  nucleation  and  growth  of  twinned  regions  are  generally  difficult  to 

predict  using  conventional  constitutive  models,  other  frameworks,  such  as  phase  field  schemes, 
have  been  used  elsewhere  (14).  The  present  FTMP-based  model  can  also  successfully  simulate 
twinning  as  demonstrated  here. 

In  the  present  simulations,  nucleation  sites  of  twinning  are  essentially  controlled  by  precursor 
inhomogeneity  due  to  slip.  More  precisely,  the  precursor  inhomogeneity  is  evolved  from  the 
FTMP  a  and  t]  terms  for  the  slip  mode,  i.e.,  Ft  (a/j ),  F2(qri)  in  equation  26,  as  shown  in  the  top 

part  of  figure  2.  This  implies  twinning  behavior  is  greatly  affected  by  the  strain  history. 
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Figure  2.  Stress-strain  responses  comparing  effect  of  precursor  inhomogeneity  on  twin  evolution. 

To  demonstrate  this  implication,  simulations  are  repeated  with  F^a^),  F2(ij/!  )  — »  0,  with  results 

shown  in  the  bottom  part  of  figure  2.  Nucleation  of  the  twin  in  this  case  is  delayed  until 
s  =  0. 163 ,  and  growth  is  ultimately  initiated  at  around  s  =  0. 170 ,  not  from  the  center  of  the 
sample  but  from  its  upper  and  lower  edges. 

Figure  3  compares  evolving  contours  of  twinned  regions  as  well  as  the  incompatibility 
term  F{Tj)Min  =  FArf  )  between  simulations  considering  the  two  conditions.  Figure  3  confirms 
that  weak  evolution  of  slip-induced  inhomogeneity,  manifested  as  retarded  growth  of  F(j] )twin , 
tends  to  inhibit  the  twin  nucleation  within  the  sample. 
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Figure  3.  Comparison  of  twin  domains  (top)  and  incompatibility  (bottom)  for  simulations  with  and  without 
dislocation-  and  incompatibility-induced  hardening. 

5.2  Flow-Evolutionary  Hypothesis  and  Duality  Diagram 

A  duality  diagram  presentation  permits  observation  of  the  associated  energy  flow  from 
elastically  stored  states  into  local  plasticity  (dissipation)  manifested  by  incompatibility  tensor 
field  evolutions.  Specifically,  such  a  diagram  is  constructed  by  plotting  the  trace  of  the 
incompatibility  tensor  tr/7  =  rjKK  versus  the  fluctuating  part  of  the  elastic  strain  energy 

8Ue  ~  SW  for  a  targeted  domain,  in  further  investigation  of  the  hypothesis  in  equation  17. 

To  focus  on  twin  nucleation  phenomena,  we  selected  an  observation  area  in  the  central  region  of 
the  domain,  such  that  the  intruding  twinning  from  the  upper  lower  edges  is  excluded,  as 
indicated  in  figure  4. 
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Figure  4.  Observation  area  for  duality  diagram  representations  based  on  flow-evolutionary  hypothesis 
comparing  cases  with  and  without  twin  nucleation. 


Figure  5  shows  duality  diagrams  for  the  observation  area  comparing  the  two  model  cases:  those 
with  and  without  twin  nucleation.  In  figure  5,  a  sharp  rise  in  the  incompatibility  ?jKK  is  observed 

as  strain  energy  Ue  increases  for  both  cases,  but  the  case  with  twin  nucleation  shows  much 
earlier  onset  of  rise  in  jjkk  than  that  without  twin  nucleation.  Examining  this  more  in  detail,  we 

see  that  the  portion  of  the  duality  diagram  where  the  sharp  rises  occur  is  enlarged  in  the  inset. 
Marked  differences  between  the  two  cases  are  clearly  confirmed.  The  case  with  twin  nucleation 
exhibits  relatively  gradual  increase  in  rjKK  with  increasing  strain  energy  Ue  before  nucleation 
takes  place  (marked  by  a  circle).  This  pregrowth  of  rjKK  corresponds  to  the  evolution  of  the 

precursor  inhomogeneity  discussed  already.  The  case  without  twin  nucleation,  in  contrast,  yields 
a  much  smaller  growth  rate  of  i]KK ,  eventually  followed  by  a  sharp  increase  not  by  nucleation 

but  by  intrusion  of  a  twinned  region  from  outside  the  observation  area. 

From  the  energetic  point  of  view,  the  duality  diagram  provides  the  following  insight.  In  the  case 
with  twin  nucleation,  a  rapid  release  of  energy  that  has  been  stored  elastically  into  the  twin 
degrees  of  freedom  results  in  twin  nucleation  and  subsequent  growth. 
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Figure  5.  Duality  diagram  for  FCC  copper  simulation  results  comparing  cases  with  and  without  twin 
nucleation  (with  and  without  dislocation-  and  incompatibility-based  hardening). 


5.3  Summary 

The  foregoing  results  demonstrate  some  prominent  capabilities  of  the  proposed  twinning  model: 

•  Natural  representations  of  twin  nucleation  are  possible  in  the  sense  that  no  artificial 
disturbance  that  promotes  the  outset  of  the  twinning  is  required. 

•  Nucleation  and  subsequent  growth  into  lenticular  shapes  is  realistically  captured. 

•  Stress-strain  responses  accompanied  by  serration  and  overall  softening  are  predicted. 

The  flow-evolutionary  examination  based  on  the  duality  diagram  representation  of  the  simulated 
results  ( tr?7  vs.  5Ue )  has  been  shown  indicative  of  twin  evolution  in  terms  of  associated  energy 
flow — specifically,  how  elastically  stored  strain  energy  is  released  to  localized  plastic  flow  or 
twinning.  This  has  been  clarified  by  comparing  simulation  results  with  and  without  the  FTMP-  ?/ 
model,  the  latter  in  which  such  twinning  evolution  does  not  occur  naturally. 
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6.  Numerical  Simulations:  HCP  Magnesium  (Mg) 


6.1  Model,  Results,  and  Discussion 

Based  on  the  kinematic  description  for  inelasticity  of  HCP  metals — three  families  of  slip  systems 
(basal,  pyramidal,  and  prismatic)  and  two  twin  systems  (tensile  and  compression  twins) — the 
FTMP-based  model  is  now  applied  to  Mg.  Basic  capabilities  for  modeling  pure  Mg,  with 
mechanical  properties  and  test  data  available  in  the  literature,  are  also  confirmed  in  this  stage, 
e.g.,  anisotropy  in  single  crystals,  focusing  on  orientation-dependent  stress-strain  responses, 
contributions  of  the  prismatic  slip  system,  tension-compression  asymmetric  stress-strain 
responses,  and  strain  rate  effects  including  rate  sensitivity  of  flow. 

Deformation  analyses  are  conducted  for  pure  single-crystal  Mg  with  the  hexagonal  crystal 
structure  (termed  here  as  hexagonal  close  packed  or  “HCP”  following  convention  even  though 
Mg  does  not  demonstrate  perfect  close  packing  [14]),  and  descriptive  capabilities  of  the  model 
are  confirmed  by  critical  comparisons  with  experimental  data  for  plane  strain  compression  in 
multiple  orientations  from  Kelley  and  Hosford  (15, 16). 

Kelley  and  Hosford  (15, 16)  conducted  a  systematic  series  of  constrained  compression 
experiments  on  Mg  single  crystals,  as  schematically  illustrated  in  figure  6.  In  the  present 
simulations,  the  die  constraint  condition  is  realistically  mimicked  by  the  plane  strain  assumption. 
Orientations  E  and  F  exhibit  peculiar  stress  responses  showing  plateau-like  stagnation  followed 
by  a  rapid  stress  increase.  The  plateau-like  stress  response  observed  for  these  two  orientations  is 
thought  to  be  an  indication  of  twin  activity  as  the  dominant  deformation  mode,  while  rapid  stress 
increase  is  only  partially  relieved  by  the  onset  of  hard  nonbasal  slip  modes  when  twinning  has 
saturated.  Orientation  G  is  particularly  soft  because  of  its  favorable  Schmid  factor  for  basal  slip, 
which  has  much  lower  strength  than  other  systems. 
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Figure  6.  Simulations  for  plane  strain  compression  of  Mg  single  crystals;  twinning  is  reported  prevalent 
in  experiments  on  crystals  in  orientations  E  and  F  (75, 16). 

Figure  7  shows  the  predicted  stress-strain  responses  for  orientations  A-G  (curves)  compared 
with  the  experimental  data  (points).  We  obtained  excellent  agreement  for  all  orientations.  In 
particular,  simulated  results  for  orientations  E  and  F  successfully  reproduce  the  rather  unusual 
and  highly  nonlinear  stress-strain  responses  induced  by  twinning. 
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Figure  7.  Compressive  stress-strain  predictions  (curves)  compared  to  experimental  data  for  Mg  (75, 16). 

6.2  Sample  Size  Dependence 

One  noteworthy  feature  of  the  present  model  is  its  ability  to  represent  absolute  scale  conveyed 
through  the  evaluation  range  of  the  second  spatial  derivative  of  strain  for  obtaining  the 
incompatibility  tensor  field  in  equation  10.  Since  the  E  and  F  orientations  exhibit  peculiar 
twin-dominated  stress-strain  responses,  they  are  expected  to  be  sensitively  affected  by  sample 
dimensions.  To  explore  this  possibility,  three  domains  are  compared:  12  x  12,  24  x  24,  and 
48  x  48  pm2.  The  associated  results  are  presented  in  figure  8,  together  with  stress-strain  curves 
comparing  the  three  sample  sizes  and  the  experimental  data  points.  What  can  be  clearly 
confirmed  is  that  the  emerging  patterns  obviously  depend  on  the  sample  size  for  the  two 
orientations.  The  stress  responses  vary  slightly  depending  on  the  sample  size,  e.g.,  the  onset  of 
the  stress  rise  for  the  E  orientation  is  controlled  by  the  sample  dimension.  The  variation  ranges 
roughly  correspond  to  those  of  the  experiments,  implying  the  scatter  of  the  experimental  data 
may  be  a  natural  consequence  of  sample  size  or  geometry. 
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Figure  8.  Stress-strain  predictions  and  experimental  data  (top)  and  simulated  strain  contours  (bottom)  for 
differently  sized  simulation  domains. 

6.3  Mesh  Size  Independence 

Another  feature  of  the  current  modeling  approach  is  its  ability  to  produce  mesh-independent 
results.  To  demonstrate  this,  we  compared  simulations  with  two  mesh  divisions,  24  x  24  and  48 
x  48  crossed  triangle  elements  with  the  same  evaluation  size  for  the  incompatibility  tensor 
calculation.  Figure  9  shows  strain  contours  for  these  domains  with  different  mesh  divisions.  The 
same  general  patterns  emerge  regardless  of  mesh  size.  Corresponding  stress-strain  curves  are 
also  displayed,  showing  almost  identical  response. 

For  comparison,  consider  some  recent  conventional  models  for  slip  and  twinning  (17, 18)  that 
appear  to  introduce  a  mode  transition  condition  between  twinning  and  plastic  slip  in  an  ad  hoc 
manner.  For  example,  one  such  model  (1  7)  assumes  a  hardening  law  as 

My)  =  /?o  for  y<  yref ,  h(y)  =  \ (y/y^  )m_1  for  y  >  y^  .  (33) 

Because  such  models  do  not  include  scale  effects  or  gradient/nonlocal  terms,  these  tend  to 
behave  similarly  regardless  of  sample  geometry  and  size,  which  can  be  physically  unrealistic. 
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Figure  9.  Stress-strain  predictions  (top)  and  strain  contours  (bottom)  for  domains  with  different  meshes. 

Phase  field  models,  on  the  other  hand,  seem  to  enable  consideration  of  size/scale  effects  since 
spatial  gradients  (e.g.,  of  order  parameter(s))  enter  the  constitutive  response.  However,  we  may 
need  to  introduce  initial  fluctuations  in  the  phase  field  (twinned  region,  in  this  case),  and, 
depending  on  the  problem,  this  step  may  affect  nucleation  sites  and  the  number  of  predicted 
equilibrium  twinned  regions.  Therefore,  this  approach  requires  a  sound  rationale. 


7.  Conclusions 


A  new  rational  model  for  deformation  twinning  has  been  developed,  extending  the  original 
FTMP  theory  of  Hasebe  for  dislocation-based  crystal  plasticity.  Twinning  degrees  of  freedom 
are  taken  into  account  via  the  FTMP-based  incompatibility  model,  which  is  further  implemented 
into  a  crystalline  plasticity  constitutive  framework  via  the  hardening  law.  The  FTMP-based 
model  is  applied  first  to  study  to  a  single  slip-oriented  FCC  crystal,  and  preliminary  simulations 
have  been  conducted  under  static  conditions  to  confirm  the  model’s  basic  capabilities. 

Simulation  results  exhibit  nucleation  and  growth  of  twinned  regions,  accompanied  by  serrated 
stress  responses  with  softening.  The  predicted  nucleation  of  twins  is  examined  in  detail  using 
duality  diagrams  in  the  context  of  Hasebe’s  flow-evolutionary  hypothesis.  Numerical  simulations 
have  also  been  performed  for  pure  single-crystal  magnesium  with  HCP  structure,  and  descriptive 
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capabilities  of  the  model  have  been  confirmed  by  critical  comparisons  with  experimental  data  for 
plane-strain  compression  of  crystals  of  multiple  orientations  available  in  the  literature.  Results 
successfully  reproduce  the  unique  stress-strain  responses  induced  by  twinning.  Activities  of  the 
various  slip  systems  and  twinning  mechanisms  for  each  orientation  have  been  predicted.  Other 
noteworthy  features  addressed  here,  but  not  captured  by  conventional  continuum  plasticity 
models,  are  sample  size  dependence  and  mesh  size  independence. 

Future  work  will  consider  dynamic  conditions  with  implementation  using  explicit  FEM. 
Hypervelocity  impact  compression  of  metals  will  be  explored,  focusing  on  substructure 
evolution,  where  experimental  observations  have  been  documented  in  the  literature. 
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